#include<iostream>
#include<vector>
#include<queue>
#include<malloc.h>
#include <math.h>
#include <lapacke.h>
#include <algorithm>
#include <iterator>


void A_M(int order)

{
  int i,j,n,q;
  double u,k1,k2,T1,T2,b4,b3,b2,b1,b0,y01,y02,y11,y12,y21,y22,y31,y32,y41,y42,y1,y2;
  double c11,c12,c13,c14,c21,c22,c23,c24,c31,c32,c33,c34,c41,c42,c43,c44,d1,d2,d3,d4;
  double *u1,*u2,*u4,*u5,v1,v2,v4,v5,tol;
  double e1,e2,emaxnorm,premaxnorm,Crate;
  double start,finish,duration;

  u1=(double *)malloc(sizeof(double)*100000000);
  u2=(double *)malloc(sizeof(double)*100000000);
  u4=(double *)malloc(sizeof(double)*100000000);
  u5=(double *)malloc(sizeof(double)*100000000);

  u1[0]=0.994;u2[0]=0;u4[0]=0;u5[0]=-2.0015851063790825224;u=1/81.45;
  T1=17.06521656015796;T2=19.14045706162071;

   if (order==1){
    b4=0.5;b3=0.5;b2=0;b1=0;b0=0;}
  else if (order==2){
    b4=5.0/12.0;b3=8.0/12.0;b2=-1.0/12.0;b1=0;b0=0;}
  else if (order==3){
    b4=9.0/24.0;b3=19.0/24.0;b2=-5.0/24.0;b1=1.0/24.0;b0=0;}
  else{
    b4=251.0/720.0;b3=646.0/720.0;b2=-264.0/720.0;b1=106.0/720.0;b0=-19.0/720.0;}

   std::cout<<"p="<<order+1<<std::endl;
      std::cout<<"the A_M method(situation1):"<<std::endl;
      std::cout<<"  节点数  "<<"  运行时间  "<<" 误差的最大范数  "<<"  收敛阶  "<<std::endl;
    
  //
 
  
  n=1e4;k1=T1/n;
  for (i=1;i<=3;i++){
    u1[i]=u1[i-1]+k1*u4[i-1];
    u2[i]=u2[i-1]+k1*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    u4[i]=u4[i-1]+k1*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    u5[i]=u5[i-1]+k1*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
  }

  
  for (i=4;i<=n;i++){
    v1=u1[i-1];
    v2=u2[i-1];
    v4=u4[i-1];
    v5=u5[i-1];

    //求右侧向量;
    d1=-k1*(b0*u4[i-4]+b1*u4[i-3]+b2*u4[i-2]+b3*u4[i-1]+b4*u4[i-1]);
    d2=-k1*(b0*u5[i-4]+b1*u5[i-3]+b2*u5[i-2]+b3*u5[i-1]+b4*u5[i-1]);
    y01=u2[i-4]*u2[i-4]+(u1[i-4]+u-1)*(u1[i-4]+u-1);y01=pow(y01,1.5);
    y11=u2[i-3]*u2[i-3]+(u1[i-3]+u-1)*(u1[i-3]+u-1);y11=pow(y11,1.5);
    y21=u2[i-2]*u2[i-2]+(u1[i-2]+u-1)*(u1[i-2]+u-1);y21=pow(y21,1.5);
    y31=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y31=pow(y31,1.5);
    y02=u2[i-4]*u2[i-4]+(u1[i-4]+u)*(u1[i-4]+u);y02=pow(y02,1.5);
    y12=u2[i-3]*u2[i-3]+(u1[i-3]+u)*(u1[i-3]+u);y12=pow(y12,1.5);
    y22=u2[i-2]*u2[i-2]+(u1[i-2]+u)*(u1[i-2]+u);y22=pow(y22,1.5);    
    y32=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y32=pow(y32,1.5);
    d3=-k1*(b0*(2*u5[i-4]+u1[i-4]-(u*(u1[i-4]+u-1))/y01-((1-u)*(u1[i-4]+u))/y02)+b1*(2*u5[i-3]+u1[i-3]-(u*(u1[i-3]+u-1))/y11-((1-u)*(u1[i-3]+u))/y12)
	   +b2*(2*u5[i-2]+u1[i-2]-(u*(u1[i-2]+u-1))/y21-((1-u)*(u1[i-2]+u))/y22)+b3*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y31-((1-u)*(u1[i-1]+u))/y32)
	    +b4*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y31-((1-u)*(u1[i-1]+u))/y32));
    d4=-k1*(b0*(u2[i-4]-2*u4[i-4]-(u*u2[i-4])/y01-((1-u)*u2[i-4])/y02)+b1*(u2[i-3]-2*u4[i-3]-(u*u2[i-3])/y11-((1-u)*u2[i-3])/y12)
	    +b2*(u2[i-2]-2*u4[i-2]-(u*u2[i-2])/y21-((1-u)*u2[i-2])/y22)+b3*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y31-((1-u)*u2[i-1])/y32)
	    +b4*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y31-((1-u)*u2[i-1])/y32));
    d1=-d1;d2=-d2;d3=-d3;d4=-d4;



    
    //求雅克比矩阵;
	    c11=1;c12=0;c13=-k1*b4;c14=0;
	    c21=0;c22=1;c23=0;c24=-k1*b4;
	    y31=pow(y31,5.0/3.0);y32=pow(y32,5.0/3.0);
	    c31=-k1*b4*(1-u*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u-1)*(u1[i-1]+u-1))/y31-(1-u)*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u)*(u1[i-1]+u))/y32);
	    c32=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y31+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y32);c33=1;c34=-2*k1*b4;
	    c41=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y31+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y32);
	    c42=-k1*b4*(1-u*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1))/y31-(1-u)*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u))/y32);
	    c43=2*k1*b4;c44=1;


	    
    //求解方程;
    double C[16]={c11,c21,c31,c41,
		  c12,c22,c32,c42,
		  c13,c23,c33,c43,
		  c14,c24,c34,c44};
    double D[4]={d1,d2,d3,d4};
    int ipiv[4];
    LAPACKE_dgesv(LAPACK_COL_MAJOR, 4, 1, C, 4, ipiv, D, 4);
    u1[i]=u1[i-1]+D[0];
    u2[i]=u2[i-1]+D[1];
    u4[i]=u4[i-1]+D[2];
    u5[i]=u5[i-1]+D[3];
    tol=fabs(D[0])+fabs(D[1])+fabs(D[2])+fabs(D[3]);
    for (j=1;j<=100;j++){
      if (tol>1e-6){
	u1[i-1]=u1[i];
	u2[i-1]=u2[i];
	u4[i-1]=u4[i];
	u5[i-1]=u5[i];
	//newton迭代
	//求右侧向量;
    d1=u1[i-1]-v1-k1*(b0*u4[i-4]+b1*u4[i-3]+b2*u4[i-2]+b3*v4+b4*u4[i-1]);
    d2=u2[i-1]-v2-k1*(b0*u5[i-4]+b1*u5[i-3]+b2*u5[i-2]+b3*v5+b4*u5[i-1]);
    y01=u2[i-4]*u2[i-4]+(u1[i-4]+u-1)*(u1[i-4]+u-1);y01=pow(y01,1.5);
    y11=u2[i-3]*u2[i-3]+(u1[i-3]+u-1)*(u1[i-3]+u-1);y11=pow(y11,1.5);
    y21=u2[i-2]*u2[i-2]+(u1[i-2]+u-1)*(u1[i-2]+u-1);y21=pow(y21,1.5);
    y31=v2*v2+(v1+u-1)*(v1+u-1);y31=pow(y31,1.5);
    y41=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y41=pow(y41,1.5);
    y02=u2[i-4]*u2[i-4]+(u1[i-4]+u)*(u1[i-4]+u);y02=pow(y02,1.5);
    y12=u2[i-3]*u2[i-3]+(u1[i-3]+u)*(u1[i-3]+u);y12=pow(y12,1.5);
    y22=u2[i-2]*u2[i-2]+(u1[i-2]+u)*(u1[i-2]+u);y22=pow(y22,1.5);    
    y32=v2*v2+(v1+u)*(v1+u);y32=pow(y32,1.5);
    y42=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y42=pow(y42,1.5);
    d3=u4[i-1]-v4-k1*(b0*(2*u5[i-4]+u1[i-4]-(u*(u1[i-4]+u-1))/y01-((1-u)*(u1[i-4]+u))/y02)+b1*(2*u5[i-3]+u1[i-3]-(u*(u1[i-3]+u-1))/y11-((1-u)*(u1[i-3]+u))/y12)
	   +b2*(2*u5[i-2]+u1[i-2]-(u*(u1[i-2]+u-1))/y21-((1-u)*(u1[i-2]+u))/y22)+b3*(2*v5+v1-(u*(v1+u-1))/y31-((1-u)*(v1+u))/y32)
	    +b4*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y41-((1-u)*(u1[i-1]+u))/y42));
    d4=u5[i-1]-v5-k1*(b0*(u2[i-4]-2*u4[i-4]-(u*u2[i-4])/y01-((1-u)*u2[i-4])/y02)+b1*(u2[i-3]-2*u4[i-3]-(u*u2[i-3])/y11-((1-u)*u2[i-3])/y12)
	    +b2*(u2[i-2]-2*u4[i-2]-(u*u2[i-2])/y21-((1-u)*u2[i-2])/y22)+b3*(v2-2*v4-(u*v2)/y31-((1-u)*v2)/y32)
		      +b4*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y41-((1-u)*u2[i-1])/y42));
    d1=-d1;d2=-d2;d3=-d3;d4=-d4;

    //求雅克比矩阵;
	    c11=1;c12=0;c13=-k1*b4;c14=0;
	    c21=0;c22=1;c23=0;c24=-k1*b4;
	    y41=pow(y41,5.0/3.0);y42=pow(y42,5.0/3.0);
	    c31=-k1*b4*(1-u*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u-1)*(u1[i-1]+u-1))/y41-(1-u)*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u)*(u1[i-1]+u))/y42);
	    c32=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y41+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y42);c33=1;c34=-2*k1*b4;
	    c41=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y41+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y42);
	    c42=-k1*b4*(1-u*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1))/y41-(1-u)*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u))/y42);
	    c43=2*k1*b4;c44=1;

    //求解方程;
    double C[16]={c11,c21,c31,c41,
		  c12,c22,c32,c42,
		  c13,c23,c33,c43,
		  c14,c24,c34,c44};
    double D[4]={d1,d2,d3,d4};
    int ipiv[4];
    LAPACKE_dgesv(LAPACK_COL_MAJOR, 4, 1, C, 4, ipiv, D, 4);
    u1[i]=u1[i-1]+D[0];
    u2[i]=u2[i-1]+D[1];
    u4[i]=u4[i-1]+D[2];
    u5[i]=u5[i-1]+D[3];
    tol=fabs(D[0])+fabs(D[1])+fabs(D[2])+fabs(D[3]);
      }
      else{
	break;}
	}
    u1[i-1]=v1;
    u2[i-1]=v2;
    u4[i-1]=v4;
    u5[i-1]=v5;
  }
  e1=u1[n]-u1[0];
  e2=u2[n]-u2[0];
  e1=fabs(e1);e2=fabs(e2);
  emaxnorm=std::max(e1,e2);

 
  //

  for (q=1;q<=10;q++){
  premaxnorm=emaxnorm;
  start = clock();    
  n=n*2;
  k1=T1/n;

  for (i=1;i<=3;i++){
    u1[i]=u1[i-1]+k1*u4[i-1];
    u2[i]=u2[i-1]+k1*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    u4[i]=u4[i-1]+k1*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    u5[i]=u5[i-1]+k1*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
  }

  
  for (i=4;i<=n;i++){
    v1=u1[i-1];
    v2=u2[i-1];
    v4=u4[i-1];
    v5=u5[i-1];

    //求右侧向量;
    d1=-k1*(b0*u4[i-4]+b1*u4[i-3]+b2*u4[i-2]+b3*u4[i-1]+b4*u4[i-1]);
    d2=-k1*(b0*u5[i-4]+b1*u5[i-3]+b2*u5[i-2]+b3*u5[i-1]+b4*u5[i-1]);
    y01=u2[i-4]*u2[i-4]+(u1[i-4]+u-1)*(u1[i-4]+u-1);y01=pow(y01,1.5);
    y11=u2[i-3]*u2[i-3]+(u1[i-3]+u-1)*(u1[i-3]+u-1);y11=pow(y11,1.5);
    y21=u2[i-2]*u2[i-2]+(u1[i-2]+u-1)*(u1[i-2]+u-1);y21=pow(y21,1.5);
    y31=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y31=pow(y31,1.5);
    y02=u2[i-4]*u2[i-4]+(u1[i-4]+u)*(u1[i-4]+u);y02=pow(y02,1.5);
    y12=u2[i-3]*u2[i-3]+(u1[i-3]+u)*(u1[i-3]+u);y12=pow(y12,1.5);
    y22=u2[i-2]*u2[i-2]+(u1[i-2]+u)*(u1[i-2]+u);y22=pow(y22,1.5);    
    y32=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y32=pow(y32,1.5);
    d3=-k1*(b0*(2*u5[i-4]+u1[i-4]-(u*(u1[i-4]+u-1))/y01-((1-u)*(u1[i-4]+u))/y02)+b1*(2*u5[i-3]+u1[i-3]-(u*(u1[i-3]+u-1))/y11-((1-u)*(u1[i-3]+u))/y12)
	   +b2*(2*u5[i-2]+u1[i-2]-(u*(u1[i-2]+u-1))/y21-((1-u)*(u1[i-2]+u))/y22)+b3*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y31-((1-u)*(u1[i-1]+u))/y32)
	    +b4*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y31-((1-u)*(u1[i-1]+u))/y32));
    d4=-k1*(b0*(u2[i-4]-2*u4[i-4]-(u*u2[i-4])/y01-((1-u)*u2[i-4])/y02)+b1*(u2[i-3]-2*u4[i-3]-(u*u2[i-3])/y11-((1-u)*u2[i-3])/y12)
	    +b2*(u2[i-2]-2*u4[i-2]-(u*u2[i-2])/y21-((1-u)*u2[i-2])/y22)+b3*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y31-((1-u)*u2[i-1])/y32)
	    +b4*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y31-((1-u)*u2[i-1])/y32));
    d1=-d1;d2=-d2;d3=-d3;d4=-d4;



    
    //求雅克比矩阵;
	    c11=1;c12=0;c13=-k1*b4;c14=0;
	    c21=0;c22=1;c23=0;c24=-k1*b4;
	    y31=pow(y31,5.0/3.0);y32=pow(y32,5.0/3.0);
	    c31=-k1*b4*(1-u*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u-1)*(u1[i-1]+u-1))/y31-(1-u)*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u)*(u1[i-1]+u))/y32);
	    c32=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y31+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y32);c33=1;c34=-2*k1*b4;
	    c41=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y31+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y32);
	    c42=-k1*b4*(1-u*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1))/y31-(1-u)*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u))/y32);
	    c43=2*k1*b4;c44=1;


	    
    //求解方程;
    double C[16]={c11,c21,c31,c41,
		  c12,c22,c32,c42,
		  c13,c23,c33,c43,
		  c14,c24,c34,c44};
    double D[4]={d1,d2,d3,d4};
    int ipiv[4];
    LAPACKE_dgesv(LAPACK_COL_MAJOR, 4, 1, C, 4, ipiv, D, 4);
    u1[i]=u1[i-1]+D[0];
    u2[i]=u2[i-1]+D[1];
    u4[i]=u4[i-1]+D[2];
    u5[i]=u5[i-1]+D[3];
    tol=fabs(D[0])+fabs(D[1])+fabs(D[2])+fabs(D[3]);
    for (j=1;j<=100;j++){
      if (tol>1e-6){
	u1[i-1]=u1[i];
	u2[i-1]=u2[i];
	u4[i-1]=u4[i];
	u5[i-1]=u5[i];
	//newton迭代
	//求右侧向量;
    d1=u1[i-1]-v1-k1*(b0*u4[i-4]+b1*u4[i-3]+b2*u4[i-2]+b3*v4+b4*u4[i-1]);
    d2=u2[i-1]-v2-k1*(b0*u5[i-4]+b1*u5[i-3]+b2*u5[i-2]+b3*v5+b4*u5[i-1]);
    y01=u2[i-4]*u2[i-4]+(u1[i-4]+u-1)*(u1[i-4]+u-1);y01=pow(y01,1.5);
    y11=u2[i-3]*u2[i-3]+(u1[i-3]+u-1)*(u1[i-3]+u-1);y11=pow(y11,1.5);
    y21=u2[i-2]*u2[i-2]+(u1[i-2]+u-1)*(u1[i-2]+u-1);y21=pow(y21,1.5);
    y31=v2*v2+(v1+u-1)*(v1+u-1);y31=pow(y31,1.5);
    y41=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y41=pow(y41,1.5);
    y02=u2[i-4]*u2[i-4]+(u1[i-4]+u)*(u1[i-4]+u);y02=pow(y02,1.5);
    y12=u2[i-3]*u2[i-3]+(u1[i-3]+u)*(u1[i-3]+u);y12=pow(y12,1.5);
    y22=u2[i-2]*u2[i-2]+(u1[i-2]+u)*(u1[i-2]+u);y22=pow(y22,1.5);    
    y32=v2*v2+(v1+u)*(v1+u);y32=pow(y32,1.5);
    y42=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y42=pow(y42,1.5);
    d3=u4[i-1]-v4-k1*(b0*(2*u5[i-4]+u1[i-4]-(u*(u1[i-4]+u-1))/y01-((1-u)*(u1[i-4]+u))/y02)+b1*(2*u5[i-3]+u1[i-3]-(u*(u1[i-3]+u-1))/y11-((1-u)*(u1[i-3]+u))/y12)
	   +b2*(2*u5[i-2]+u1[i-2]-(u*(u1[i-2]+u-1))/y21-((1-u)*(u1[i-2]+u))/y22)+b3*(2*v5+v1-(u*(v1+u-1))/y31-((1-u)*(v1+u))/y32)
	    +b4*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y41-((1-u)*(u1[i-1]+u))/y42));
    d4=u5[i-1]-v5-k1*(b0*(u2[i-4]-2*u4[i-4]-(u*u2[i-4])/y01-((1-u)*u2[i-4])/y02)+b1*(u2[i-3]-2*u4[i-3]-(u*u2[i-3])/y11-((1-u)*u2[i-3])/y12)
	    +b2*(u2[i-2]-2*u4[i-2]-(u*u2[i-2])/y21-((1-u)*u2[i-2])/y22)+b3*(v2-2*v4-(u*v2)/y31-((1-u)*v2)/y32)
		      +b4*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y41-((1-u)*u2[i-1])/y42));
    d1=-d1;d2=-d2;d3=-d3;d4=-d4;

    //求雅克比矩阵;
	    c11=1;c12=0;c13=-k1*b4;c14=0;
	    c21=0;c22=1;c23=0;c24=-k1*b4;
	    y41=pow(y41,5.0/3.0);y42=pow(y42,5.0/3.0);
	    c31=-k1*b4*(1-u*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u-1)*(u1[i-1]+u-1))/y41-(1-u)*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u)*(u1[i-1]+u))/y42);
	    c32=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y41+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y42);c33=1;c34=-2*k1*b4;
	    c41=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y41+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y42);
	    c42=-k1*b4*(1-u*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1))/y41-(1-u)*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u))/y42);
	    c43=2*k1*b4;c44=1;

    //求解方程;
    double C[16]={c11,c21,c31,c41,
		  c12,c22,c32,c42,
		  c13,c23,c33,c43,
		  c14,c24,c34,c44};
    double D[4]={d1,d2,d3,d4};
    int ipiv[4];
    LAPACKE_dgesv(LAPACK_COL_MAJOR, 4, 1, C, 4, ipiv, D, 4);
    u1[i]=u1[i-1]+D[0];
    u2[i]=u2[i-1]+D[1];
    u4[i]=u4[i-1]+D[2];
    u5[i]=u5[i-1]+D[3];
    tol=fabs(D[0])+fabs(D[1])+fabs(D[2])+fabs(D[3]);
      }
      else{
	break;}
	}
    u1[i-1]=v1;
    u2[i-1]=v2;
    u4[i-1]=v4;
    u5[i-1]=v5;
  }
  e1=u1[n]-u1[0];
  e2=u2[n]-u2[0];
  e1=fabs(e1);e2=fabs(e2);
  emaxnorm=std::max(e1,e2);
  Crate=premaxnorm/emaxnorm;Crate=log2(Crate);
  finish = clock();
  duration = (double)(finish - start) / CLOCKS_PER_SEC;
  std::cout<<"  "<<n<<"     " <<duration<<"    "  <<emaxnorm<<"      "<<Crate  <<std::endl;
 
  }

  T1=T2,u1[0]=0.87978;u2[0]=0;u4[0]=0;u5[0]=-0.3797;
      std::cout<<"the A_M method(situation2):"<<std::endl;
      std::cout<<"  节点数  "<<"  运行时间  "<<" 误差的最大范数  "<<"  收敛阶  "<<std::endl;

  n=5e2;k1=T1/n;
  for (i=1;i<=3;i++){
    u1[i]=u1[i-1]+k1*u4[i-1];
    u2[i]=u2[i-1]+k1*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    u4[i]=u4[i-1]+k1*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    u5[i]=u5[i-1]+k1*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
  }

  
  for (i=4;i<=n;i++){
    v1=u1[i-1];
    v2=u2[i-1];
    v4=u4[i-1];
    v5=u5[i-1];

    //求右侧向量;
    d1=-k1*(b0*u4[i-4]+b1*u4[i-3]+b2*u4[i-2]+b3*u4[i-1]+b4*u4[i-1]);
    d2=-k1*(b0*u5[i-4]+b1*u5[i-3]+b2*u5[i-2]+b3*u5[i-1]+b4*u5[i-1]);
    y01=u2[i-4]*u2[i-4]+(u1[i-4]+u-1)*(u1[i-4]+u-1);y01=pow(y01,1.5);
    y11=u2[i-3]*u2[i-3]+(u1[i-3]+u-1)*(u1[i-3]+u-1);y11=pow(y11,1.5);
    y21=u2[i-2]*u2[i-2]+(u1[i-2]+u-1)*(u1[i-2]+u-1);y21=pow(y21,1.5);
    y31=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y31=pow(y31,1.5);
    y02=u2[i-4]*u2[i-4]+(u1[i-4]+u)*(u1[i-4]+u);y02=pow(y02,1.5);
    y12=u2[i-3]*u2[i-3]+(u1[i-3]+u)*(u1[i-3]+u);y12=pow(y12,1.5);
    y22=u2[i-2]*u2[i-2]+(u1[i-2]+u)*(u1[i-2]+u);y22=pow(y22,1.5);    
    y32=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y32=pow(y32,1.5);
    d3=-k1*(b0*(2*u5[i-4]+u1[i-4]-(u*(u1[i-4]+u-1))/y01-((1-u)*(u1[i-4]+u))/y02)+b1*(2*u5[i-3]+u1[i-3]-(u*(u1[i-3]+u-1))/y11-((1-u)*(u1[i-3]+u))/y12)
	   +b2*(2*u5[i-2]+u1[i-2]-(u*(u1[i-2]+u-1))/y21-((1-u)*(u1[i-2]+u))/y22)+b3*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y31-((1-u)*(u1[i-1]+u))/y32)
	    +b4*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y31-((1-u)*(u1[i-1]+u))/y32));
    d4=-k1*(b0*(u2[i-4]-2*u4[i-4]-(u*u2[i-4])/y01-((1-u)*u2[i-4])/y02)+b1*(u2[i-3]-2*u4[i-3]-(u*u2[i-3])/y11-((1-u)*u2[i-3])/y12)
	    +b2*(u2[i-2]-2*u4[i-2]-(u*u2[i-2])/y21-((1-u)*u2[i-2])/y22)+b3*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y31-((1-u)*u2[i-1])/y32)
	    +b4*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y31-((1-u)*u2[i-1])/y32));
    d1=-d1;d2=-d2;d3=-d3;d4=-d4;



    
    //求雅克比矩阵;
	    c11=1;c12=0;c13=-k1*b4;c14=0;
	    c21=0;c22=1;c23=0;c24=-k1*b4;
	    y31=pow(y31,5.0/3.0);y32=pow(y32,5.0/3.0);
	    c31=-k1*b4*(1-u*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u-1)*(u1[i-1]+u-1))/y31-(1-u)*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u)*(u1[i-1]+u))/y32);
	    c32=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y31+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y32);c33=1;c34=-2*k1*b4;
	    c41=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y31+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y32);
	    c42=-k1*b4*(1-u*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1))/y31-(1-u)*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u))/y32);
	    c43=2*k1*b4;c44=1;


	    
    //求解方程;
    double C[16]={c11,c21,c31,c41,
		  c12,c22,c32,c42,
		  c13,c23,c33,c43,
		  c14,c24,c34,c44};
    double D[4]={d1,d2,d3,d4};
    int ipiv[4];
    LAPACKE_dgesv(LAPACK_COL_MAJOR, 4, 1, C, 4, ipiv, D, 4);
    u1[i]=u1[i-1]+D[0];
    u2[i]=u2[i-1]+D[1];
    u4[i]=u4[i-1]+D[2];
    u5[i]=u5[i-1]+D[3];
    tol=fabs(D[0])+fabs(D[1])+fabs(D[2])+fabs(D[3]);
    for (j=1;j<=100;j++){
      if (tol>1e-6){
	u1[i-1]=u1[i];
	u2[i-1]=u2[i];
	u4[i-1]=u4[i];
	u5[i-1]=u5[i];
	//newton迭代
	//求右侧向量;
    d1=u1[i-1]-v1-k1*(b0*u4[i-4]+b1*u4[i-3]+b2*u4[i-2]+b3*v4+b4*u4[i-1]);
    d2=u2[i-1]-v2-k1*(b0*u5[i-4]+b1*u5[i-3]+b2*u5[i-2]+b3*v5+b4*u5[i-1]);
    y01=u2[i-4]*u2[i-4]+(u1[i-4]+u-1)*(u1[i-4]+u-1);y01=pow(y01,1.5);
    y11=u2[i-3]*u2[i-3]+(u1[i-3]+u-1)*(u1[i-3]+u-1);y11=pow(y11,1.5);
    y21=u2[i-2]*u2[i-2]+(u1[i-2]+u-1)*(u1[i-2]+u-1);y21=pow(y21,1.5);
    y31=v2*v2+(v1+u-1)*(v1+u-1);y31=pow(y31,1.5);
    y41=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y41=pow(y41,1.5);
    y02=u2[i-4]*u2[i-4]+(u1[i-4]+u)*(u1[i-4]+u);y02=pow(y02,1.5);
    y12=u2[i-3]*u2[i-3]+(u1[i-3]+u)*(u1[i-3]+u);y12=pow(y12,1.5);
    y22=u2[i-2]*u2[i-2]+(u1[i-2]+u)*(u1[i-2]+u);y22=pow(y22,1.5);    
    y32=v2*v2+(v1+u)*(v1+u);y32=pow(y32,1.5);
    y42=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y42=pow(y42,1.5);
    d3=u4[i-1]-v4-k1*(b0*(2*u5[i-4]+u1[i-4]-(u*(u1[i-4]+u-1))/y01-((1-u)*(u1[i-4]+u))/y02)+b1*(2*u5[i-3]+u1[i-3]-(u*(u1[i-3]+u-1))/y11-((1-u)*(u1[i-3]+u))/y12)
	   +b2*(2*u5[i-2]+u1[i-2]-(u*(u1[i-2]+u-1))/y21-((1-u)*(u1[i-2]+u))/y22)+b3*(2*v5+v1-(u*(v1+u-1))/y31-((1-u)*(v1+u))/y32)
	    +b4*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y41-((1-u)*(u1[i-1]+u))/y42));
    d4=u5[i-1]-v5-k1*(b0*(u2[i-4]-2*u4[i-4]-(u*u2[i-4])/y01-((1-u)*u2[i-4])/y02)+b1*(u2[i-3]-2*u4[i-3]-(u*u2[i-3])/y11-((1-u)*u2[i-3])/y12)
	    +b2*(u2[i-2]-2*u4[i-2]-(u*u2[i-2])/y21-((1-u)*u2[i-2])/y22)+b3*(v2-2*v4-(u*v2)/y31-((1-u)*v2)/y32)
		      +b4*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y41-((1-u)*u2[i-1])/y42));
    d1=-d1;d2=-d2;d3=-d3;d4=-d4;

    //求雅克比矩阵;
	    c11=1;c12=0;c13=-k1*b4;c14=0;
	    c21=0;c22=1;c23=0;c24=-k1*b4;
	    y41=pow(y41,5.0/3.0);y42=pow(y42,5.0/3.0);
	    c31=-k1*b4*(1-u*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u-1)*(u1[i-1]+u-1))/y41-(1-u)*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u)*(u1[i-1]+u))/y42);
	    c32=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y41+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y42);c33=1;c34=-2*k1*b4;
	    c41=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y41+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y42);
	    c42=-k1*b4*(1-u*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1))/y41-(1-u)*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u))/y42);
	    c43=2*k1*b4;c44=1;

    //求解方程;
    double C[16]={c11,c21,c31,c41,
		  c12,c22,c32,c42,
		  c13,c23,c33,c43,
		  c14,c24,c34,c44};
    double D[4]={d1,d2,d3,d4};
    int ipiv[4];
    LAPACKE_dgesv(LAPACK_COL_MAJOR, 4, 1, C, 4, ipiv, D, 4);
    u1[i]=u1[i-1]+D[0];
    u2[i]=u2[i-1]+D[1];
    u4[i]=u4[i-1]+D[2];
    u5[i]=u5[i-1]+D[3];
    tol=fabs(D[0])+fabs(D[1])+fabs(D[2])+fabs(D[3]);
      }
      else{
	break;}
	}
    u1[i-1]=v1;
    u2[i-1]=v2;
    u4[i-1]=v4;
    u5[i-1]=v5;
  }
  e1=u1[n]-u1[0];
  e2=u2[n]-u2[0];
  e1=fabs(e1);e2=fabs(e2);
  emaxnorm=std::max(e1,e2);

 
  //

  for (q=1;q<=10;q++){
  premaxnorm=emaxnorm;
  start = clock();    
  n=n*2;
  k1=T1/n;

  for (i=1;i<=3;i++){
    u1[i]=u1[i-1]+k1*u4[i-1];
    u2[i]=u2[i-1]+k1*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    u4[i]=u4[i-1]+k1*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    u5[i]=u5[i-1]+k1*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
  }

  
  for (i=4;i<=n;i++){
    v1=u1[i-1];
    v2=u2[i-1];
    v4=u4[i-1];
    v5=u5[i-1];

    //求右侧向量;
    d1=-k1*(b0*u4[i-4]+b1*u4[i-3]+b2*u4[i-2]+b3*u4[i-1]+b4*u4[i-1]);
    d2=-k1*(b0*u5[i-4]+b1*u5[i-3]+b2*u5[i-2]+b3*u5[i-1]+b4*u5[i-1]);
    y01=u2[i-4]*u2[i-4]+(u1[i-4]+u-1)*(u1[i-4]+u-1);y01=pow(y01,1.5);
    y11=u2[i-3]*u2[i-3]+(u1[i-3]+u-1)*(u1[i-3]+u-1);y11=pow(y11,1.5);
    y21=u2[i-2]*u2[i-2]+(u1[i-2]+u-1)*(u1[i-2]+u-1);y21=pow(y21,1.5);
    y31=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y31=pow(y31,1.5);
    y02=u2[i-4]*u2[i-4]+(u1[i-4]+u)*(u1[i-4]+u);y02=pow(y02,1.5);
    y12=u2[i-3]*u2[i-3]+(u1[i-3]+u)*(u1[i-3]+u);y12=pow(y12,1.5);
    y22=u2[i-2]*u2[i-2]+(u1[i-2]+u)*(u1[i-2]+u);y22=pow(y22,1.5);    
    y32=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y32=pow(y32,1.5);
    d3=-k1*(b0*(2*u5[i-4]+u1[i-4]-(u*(u1[i-4]+u-1))/y01-((1-u)*(u1[i-4]+u))/y02)+b1*(2*u5[i-3]+u1[i-3]-(u*(u1[i-3]+u-1))/y11-((1-u)*(u1[i-3]+u))/y12)
	   +b2*(2*u5[i-2]+u1[i-2]-(u*(u1[i-2]+u-1))/y21-((1-u)*(u1[i-2]+u))/y22)+b3*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y31-((1-u)*(u1[i-1]+u))/y32)
	    +b4*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y31-((1-u)*(u1[i-1]+u))/y32));
    d4=-k1*(b0*(u2[i-4]-2*u4[i-4]-(u*u2[i-4])/y01-((1-u)*u2[i-4])/y02)+b1*(u2[i-3]-2*u4[i-3]-(u*u2[i-3])/y11-((1-u)*u2[i-3])/y12)
	    +b2*(u2[i-2]-2*u4[i-2]-(u*u2[i-2])/y21-((1-u)*u2[i-2])/y22)+b3*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y31-((1-u)*u2[i-1])/y32)
	    +b4*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y31-((1-u)*u2[i-1])/y32));
    d1=-d1;d2=-d2;d3=-d3;d4=-d4;



    
    //求雅克比矩阵;
	    c11=1;c12=0;c13=-k1*b4;c14=0;
	    c21=0;c22=1;c23=0;c24=-k1*b4;
	    y31=pow(y31,5.0/3.0);y32=pow(y32,5.0/3.0);
	    c31=-k1*b4*(1-u*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u-1)*(u1[i-1]+u-1))/y31-(1-u)*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u)*(u1[i-1]+u))/y32);
	    c32=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y31+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y32);c33=1;c34=-2*k1*b4;
	    c41=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y31+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y32);
	    c42=-k1*b4*(1-u*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1))/y31-(1-u)*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u))/y32);
	    c43=2*k1*b4;c44=1;


	    
    //求解方程;
    double C[16]={c11,c21,c31,c41,
		  c12,c22,c32,c42,
		  c13,c23,c33,c43,
		  c14,c24,c34,c44};
    double D[4]={d1,d2,d3,d4};
    int ipiv[4];
    LAPACKE_dgesv(LAPACK_COL_MAJOR, 4, 1, C, 4, ipiv, D, 4);
    u1[i]=u1[i-1]+D[0];
    u2[i]=u2[i-1]+D[1];
    u4[i]=u4[i-1]+D[2];
    u5[i]=u5[i-1]+D[3];
    tol=fabs(D[0])+fabs(D[1])+fabs(D[2])+fabs(D[3]);
    for (j=1;j<=100;j++){
      if (tol>1e-6){
	u1[i-1]=u1[i];
	u2[i-1]=u2[i];
	u4[i-1]=u4[i];
	u5[i-1]=u5[i];
	//newton迭代
	//求右侧向量;
    d1=u1[i-1]-v1-k1*(b0*u4[i-4]+b1*u4[i-3]+b2*u4[i-2]+b3*v4+b4*u4[i-1]);
    d2=u2[i-1]-v2-k1*(b0*u5[i-4]+b1*u5[i-3]+b2*u5[i-2]+b3*v5+b4*u5[i-1]);
    y01=u2[i-4]*u2[i-4]+(u1[i-4]+u-1)*(u1[i-4]+u-1);y01=pow(y01,1.5);
    y11=u2[i-3]*u2[i-3]+(u1[i-3]+u-1)*(u1[i-3]+u-1);y11=pow(y11,1.5);
    y21=u2[i-2]*u2[i-2]+(u1[i-2]+u-1)*(u1[i-2]+u-1);y21=pow(y21,1.5);
    y31=v2*v2+(v1+u-1)*(v1+u-1);y31=pow(y31,1.5);
    y41=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y41=pow(y41,1.5);
    y02=u2[i-4]*u2[i-4]+(u1[i-4]+u)*(u1[i-4]+u);y02=pow(y02,1.5);
    y12=u2[i-3]*u2[i-3]+(u1[i-3]+u)*(u1[i-3]+u);y12=pow(y12,1.5);
    y22=u2[i-2]*u2[i-2]+(u1[i-2]+u)*(u1[i-2]+u);y22=pow(y22,1.5);    
    y32=v2*v2+(v1+u)*(v1+u);y32=pow(y32,1.5);
    y42=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y42=pow(y42,1.5);
    d3=u4[i-1]-v4-k1*(b0*(2*u5[i-4]+u1[i-4]-(u*(u1[i-4]+u-1))/y01-((1-u)*(u1[i-4]+u))/y02)+b1*(2*u5[i-3]+u1[i-3]-(u*(u1[i-3]+u-1))/y11-((1-u)*(u1[i-3]+u))/y12)
	   +b2*(2*u5[i-2]+u1[i-2]-(u*(u1[i-2]+u-1))/y21-((1-u)*(u1[i-2]+u))/y22)+b3*(2*v5+v1-(u*(v1+u-1))/y31-((1-u)*(v1+u))/y32)
	    +b4*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y41-((1-u)*(u1[i-1]+u))/y42));
    d4=u5[i-1]-v5-k1*(b0*(u2[i-4]-2*u4[i-4]-(u*u2[i-4])/y01-((1-u)*u2[i-4])/y02)+b1*(u2[i-3]-2*u4[i-3]-(u*u2[i-3])/y11-((1-u)*u2[i-3])/y12)
	    +b2*(u2[i-2]-2*u4[i-2]-(u*u2[i-2])/y21-((1-u)*u2[i-2])/y22)+b3*(v2-2*v4-(u*v2)/y31-((1-u)*v2)/y32)
		      +b4*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y41-((1-u)*u2[i-1])/y42));
    d1=-d1;d2=-d2;d3=-d3;d4=-d4;

    //求雅克比矩阵;
	    c11=1;c12=0;c13=-k1*b4;c14=0;
	    c21=0;c22=1;c23=0;c24=-k1*b4;
	    y41=pow(y41,5.0/3.0);y42=pow(y42,5.0/3.0);
	    c31=-k1*b4*(1-u*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u-1)*(u1[i-1]+u-1))/y41-(1-u)*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u)*(u1[i-1]+u))/y42);
	    c32=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y41+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y42);c33=1;c34=-2*k1*b4;
	    c41=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y41+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y42);
	    c42=-k1*b4*(1-u*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1))/y41-(1-u)*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u))/y42);
	    c43=2*k1*b4;c44=1;

    //求解方程;
    double C[16]={c11,c21,c31,c41,
		  c12,c22,c32,c42,
		  c13,c23,c33,c43,
		  c14,c24,c34,c44};
    double D[4]={d1,d2,d3,d4};
    int ipiv[4];
    LAPACKE_dgesv(LAPACK_COL_MAJOR, 4, 1, C, 4, ipiv, D, 4);
    u1[i]=u1[i-1]+D[0];
    u2[i]=u2[i-1]+D[1];
    u4[i]=u4[i-1]+D[2];
    u5[i]=u5[i-1]+D[3];
    tol=fabs(D[0])+fabs(D[1])+fabs(D[2])+fabs(D[3]);
      }
      else{
	break;}
	}
    u1[i-1]=v1;
    u2[i-1]=v2;
    u4[i-1]=v4;
    u5[i-1]=v5;
  }
  e1=u1[n]-u1[0];
  e2=u2[n]-u2[0];
  e1=fabs(e1);e2=fabs(e2);
  emaxnorm=std::max(e1,e2);
  Crate=premaxnorm/emaxnorm;Crate=log2(Crate);
  finish = clock();
  duration = (double)(finish - start) / CLOCKS_PER_SEC;
  std::cout<<"  "<<n<<"     " <<duration<<"    "  <<emaxnorm<<"      "<<Crate  <<std::endl;
 
  }

      
}

